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Abstract 
l/"~) ' We present the first example of a phase transition in a nonequilibrium steady-state that can 

o 



be argued analytically to be first order. The system of interest is a two-species reaction- 
diffusion problem whose control parameter is the total density p. Mean-field theory pre- 
0^ \ diets a second-order transition between two stationary states at a critical density p = p c . 

We develop a phenomenological picture that, instead, below the upper critical dimension 
d c = 4, predicts a first-order transition. This picture is confirmed by hysteresis found in 
numerical simulations, and by the study of a renormalization-group improved equation of 
state. The latter approach is inspired by the Weinberg-Coleman mechanism in QED. 
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1 Motivation 

1.1 Fluctuation induced first order transitions 

In the realm of equilibrium critical phenomena it is well-known that systems which 
in high space dimension d undergo a first-order transition, may exhibit a second- 
order transition below their upper critical dimension d = d c . Examples are spin 
systems with cubic anisotropy [111, EL RI], type-II superconductors [01], and the three- 
state Potts model [|5|]. In equilibrium phenomena the phase diagram can be deduced 
from the analysis of the global extrema of a free energy functional. The global free 
energy minima correspond to stable phases. For d < d c this functional has to 
incorporate fluctuation effects. When fluctuations modify the energy landscape 
to the point of changing a second-order transition into a first-order one, one has a 
fluctuation-induced first-order transition. This phenomenon is also said to be due to 
the Coleman-Weinberg mechanism [||]. Indeed formally similar phenomena were 
first found in the study of the coupling of QED radiative corrections to a charged 
scalar field. 

In nonequilibrium systems there is no such concept as a free energy; the steady 
state phase diagram cannot be deduced from a thermodynamic potential, and its 
statistical mechanics is not based on a partition function. The starting point, in- 
stead, is usually an evolution equation (often a master equation), whose stationary 
solutions are to be determined and yield the steady-state phase diagram. This indi- 
rect definition of the phase diagram renders analytic approaches very cumbersome. 
In the past twenty years techniques have been devised to find the steady states of 
such master equations and extract from them physical properties of interest (order 
parameter, correlation functions,. . . ). These techniques include short time series 
expansion, numerical simulations, real-space renormalization and field-theoretic 
approaches. Nonequilibrium steady states (NESS) may undergo phase transitions 
in the same way as do equilibrium states. Several examples of continuous transi- 
tions in such systems have been found and well studied. To our knowledge, the 
only first-order transitions in NESS known today occur in the asymmetric exclu- 
sion model in one space dimension, as demonstrated analytically in Ref. [Q]. This 
model belongs to the subclass of driven diffusive systems: due to an externally ap- 
plied field these systems have a spatially anisotropic current carrying NESS. The 
occurrence of first order transitions in certain other driven diffusive systems [|| is 
also suggested by numerical simulations. Finally, Schmittmann and Janssen [Q] 
have argued field-theoretically that a similar fluctuation mechanism may induce a 
first-order longitudinal transition (high and low density stripes perpendicular to the 
driving field) in a driven diffusive system with a single conserved density. 



The present study bears on the phase transition in the NESS of a diffusion- 
limited reaction between two species of particles. This system does not belong 
to the subclass of driven systems: it remains spatially isotropic under all circum- 
stances. 

We exhibit here, for the first time with analytic arguments, a fluctuation-induced 
first order transition in the steady-state of a reaction-diffusion process. We present 
an analytic procedure that allows to access the phase diagram of the system in a 
very explicit fashion. The outline of the article is as follows. In the next subsection 
we define the two-species reaction-diffusion model. In section 2 we recall some 
known properties of its phase diagram. We also introduce the field-theoretical for- 
malism, which will be our main tool of analysis. Section 3 presents a heuristic 
approach to the first-order transition in terms of a nucleating and diffusing droplet 
picture. In section 4 we report on numerical simulations on a two-dimensional 
system, in which a hysteresis loop is observed for the order parameter. This con- 
firms the suspected existence of a first order transition below the critical dimension 
d c = 4. Sections 5 through 8 contain the field-theoretic approach: the derivation of 
a renormalization group improved equation of state, valid in dimension d = d c — e, 
followed by the study of its solutions and of their stability with respect to spatial 
perturbations. We conclude with a series of possible applications. 

1.2 Reaction-diffusion model 

Particles of two species, A and B, diffuse in a ti-dimensional space with diffusion 
constants Da and Db, respectively. Upon encounter an A and a B are converted 
into two -B's at a rate k$ per unit of volume, 

A + B H B + B (1) 

Besides a B spontaneously decays into an A at a rate 7, 

B -^ A (2) 

Denoting the local A and B densities by pa and ps, respectively, we can write the 
mean-field equations as 

d t pB = D B Ap B - JPB + hPAPB 

d t p A = D A ApA + 1PB ~ koPAPB (3) 

The total particle density p, is a conserved quantity and will be the control param- 
eter. In the initial state particles are distributed randomly and independently, with 



a given fraction of each species. Let p^ and p^ be the steady state values of the A 
and B particle density, respectively. Obviously their sum is equal to p. One easily 
derives from (Q) that there exists a threshold density p c = -y/k such that for p > p c 
the steady state of the system is active, that is, has p B > 0, and that for p < p c it 
is absorbing, that is, has p s g < 0. Hence p s g is the order parameter for this system, 
and an important question is how this quantity behaves at the transition point. 

This model may be cast into the form of a field theory (and then turns out to 
generalize the field theory of the Directed Percolation problem). It was shown by 



field-theoretic methods in [10] and [11] that for < Da < Db the transition 



between the steady states at p = p c is continuous. It is characterized by a set of 
critical exponents that differ from their mean-field values in d < d c = 4. Whereas 
for Da < Db the phase diagram may be obtained by convential tools of analysis, 
such as renormalization group approaches based on a field-theoretic formulation of 
the dynamics, the case Da > Db cannot be analyzed along the same lines. In tech- 
nical terms, the renormalization group flows away to a region where the theory is 
ill-defined. We inteipret this as meaning that there is no continuous transition, and 
a natural idea, if one still believes in the existence of a transition, is the occurrence 
of a first-order one. The article is concerned with the case Da > Db- 

2 Field-theoretic formulation 

2.1 Langevin equations 

There are at least two different ways to construct a field theory that describes a 
reaction-diffusion problem such as the one we just defined. One of them is to en- 
code the stochastic rules for particle diffusion and reaction in a master equation 
for the probability of occurrence of a state of given local particle numbers at a 
given time. The master equation may be converted into an exactly equivalent field 



theory following methods that were pioneered by Peliti [ ]12[ ] and others. We will 



follow a different way of proceeding that has been widely employed, in particular, 



by Janssen and co-workers Q13Q. We postulate for the space and time dependent 
densities of the A and B particles two Langevin equations in which the noise terms 
have been deduced by heuristic considerations. The result of this approach differs 
from Peliti's up to terms that are irrelevant in the limit of large times and distances. 

We switch now to the notation of field theory and denote by ^(r, t) the coarse- 
grained B particle density and by m(r, t) the deviation from average of the coarse- 
grained total particle density. The deterministic part of the Langevin equations to 
be constructed should be the conventional mean-field reaction-diffusion PDE's of 
equation (||). Upon adding two noise terms r] and £ we get in the new notation, and 



after redefinition of several parameters, 



*S_,2 



dttp = XAip — Xrip ip ~~ Xfrmp + t] (4) 



d t m = Am + XaAip + £ (5) 

Here A = Db/Da is the ratio of the diffusion constants; r is proportional to the 
deviation of the total density from its mean-field critical value 7/fe; g, g and / are 
proportional to the contamination rate k; and the parameter a, which will play a 
key role in the present study, is proportional to 1 — jf-. 

For mathematical convenience we will want 77 and ^ to represent mutually un- 
corrected Gaussian white noise. The noise rj in the equation for tp should vanish 
with ip since ip = is an absorbing state. For the autocorrelation of r\ we therefore 
retain the first term of a hypothetical series expansion in powers of ip. This proce- 



dure is best described in [ |13| , |10| ]. The autocorrelation of £ should be such that m 
is locally conserved. With these conditions the simplest possible expressions for 
the autocorrelations of rj and £ are, explicitly, 

( V (r, tW, 0) = A^(r, t)SF> (r - r')8(t - t') 
(£(r, tW, t')} = 2V r V r ,5^ (r - r')S(t - t') (6) 

The Langevin equation (Q) is to be understood with the Ito discretization rule. It is 
also possible to derive these equations ab initio by the operator formalism used in 



[11]. 



Using the Janssen-De Dominicis formalism JT4] , |I3| ] we obtain the physical 
observables as functional integrals over four fields ip,Tp,m,m weighted by a factor 

exp(— S[tp, ip, m, m\), with 

S[ip, tp, m,m}= d d xdt tp(d t + A(r - A))ip + m(d t - A)m 

J - (7) 

— (Vm) 2 — AcrmA^ -\ ip 2 ^ iptfj 2 + Xfipipm 



It is possible (and sometimes more practical) to eliminate the fluctuating density 
field m and its response field m, which yields an effective action for the -0, ip fields 
alone. 



2.2 Mean-field equation of state 

Our starting point is the action describing the dynamics of the system in the pres- 
ence of an arbitrary source of B particles A/i(r, t) in which the fluctuating density 
m and its response field m have been integrated out. It reads 

S$,ip] = i d d rdt \i>{d t + \{t - A))V> + ^Vv 5 - ^W 2 - A/m/S 



2 

d d rd d r'dtdt' 



££-W(r, t)C (r -r'-t- t')^(r', t>) (8) 

- A 2 a/^(r, *)G (r - r'; i - t / )A^(r / , i') 

where the spatial Fourier transforms of Go and Co are 

G (q; t-t') = 0(t - t / )e-« a (*- t ') J C (q; t - i') = e^''^ (9) 

We first look for an a priori inhomogeneous steady state (in terms of Fourier trans- 
forms, one takes the limit u> —> 0) then one specializes the study to a homogeneous 
steady state (and one takes the limit q —> 0). The limit of infinite times (cor- 
responding to a system reaching a steady-state) and the limit of a homogeneous 
system do not commute. In the limit of a vanishing source term the mean-field 
equation of state for a homogeneous order parameter $ is found by imposing that 

lim lim -=- -[0, *] = A* ( r + ftf) = (10) 

q^o^o^(q, w ) 1 V 2 / 

with g = g — 2\o~f. It is important to note that written in terms of the original 
parameters g > as long as j^- > (the cases Db = or Da = would 
require a separate study). One may conclude that the steady state is active (a finite 
fraction of B's survive indefinitely) if r > 0, while the system eventually falls into 
an absorbing i?-free state if r < 0. Mean-field therefore predicts a continuous 
transition between those states at r = 0, independently of the ratio of the diffusion 
constants. 

2.3 Renormalization 

In order to go beyond mean-field we have performed a one-loop perturbation ex- 
pansion of the two and three-leg vertex functions. Renormalization is then required 
to extract physically relevant information from this expansion. We shall proceed 
within the framework of dimensional regularization and of the minimal subtraction 



scheme. In order to absorb the e-poles in the vertex functions into a reparametriza- 
tion of coupling constants and fields we introduce the renormalized quantities ip R , 
A R , p R etc. defined by 

j = z 1 J 2 $ R ^ = 4 /2 Vr z = (z^) 1 ' 2 

ZX = Z\\ A/9 = A R p R Zj \a = A r o-r 
Z x t = Z t t Al/ 2 Z x f = Z f U^I 2 Z x h = Z l J 2 h R (11) 
A l J 2 Z x z]l 2 g = <J R p £ ' 2 {Z g g R + Wf R ) A 1 /' 'Z x zf ''a R ~g = Zfg R p £ / 2 . 

Here p denotes an external momentum scale. The renormalization factors depend 
on u = g R g R , v = f 2 and w = f R g R and are at one-loop order given by 



u 2v 3 + p 

~4~s ~ e{l+p) 2 ~ 2e{l + pf 



Z = l + ^-^r-^-^rr-^r 2W (12) 



v 1 , u 2 P V P 2 + 4 P -1 

Zx = 1 + 8S-7(TTW- 4 £ (i + P )3 w (13) 



u 2(3 + /)) 2(2 + p) 

e e(l + p) 2V e(l + P f 
u 2(3 + p) 5 + 3p 

1 ~ e(l + p) 2V ~ e(l + P f 



Z '9 = l + -- 7JT-^2 V ~ 77T-^2 W < 14 > 



Z 9 = l + -~ 7IT-^2 V ~ 7<T-^2 W W 



Zf = \-\ ; -77 ; r^W (16) 

' 2e e{l + p) 2 e{l + p) 2 

\v 2w 

W= y—^ + j— ~ v (17) 

£p(l + p) ep(l + p) 

Since only Z is fixed by the renormalization conditions but not the individual fac- 
tors Zj, and Z^p we may set Z^ = 1. 



2.4 Renormalization group and fixed points 

^From the above Z-factors and the definition of the renormalized couplings one 
finds the flow equations for the renormalized couplings. These read 



Pu = P 



du 
dp 



u 



bare 



-e + 



3u 2(5 + 5p + 2p 2 )v 2(4 + 5p + 2p 2 )w 



a+py 



+w 



(1 + p) 3 

4v 2w 

+ 



a dv 



v(-e + 2k) 



bare 



2>U 



p{i + p) p{i + p) 

Av (9 + 8p + 3p 2 )w 



Pw= P -7- 

dp 



w 



bare 



9u 



1 + p) 3 2(l + p)3 

2(A + 2p + p 2 )v 3(9 + 8p + 3p 2 )w 



P = P 



dp 

d// 



-Cp 



(i + p) 3 

u 2v 



bare 



4(1 + pY 
(7 + 4p + p 2 )w 



and the Wilson function is 

dlnZ 



7 = P 



dp 



2v 



(1+p) 3 4(l+p)3 

(3 + p)w 



bare 



U 

"4 + (l + p) 2 ' 2(1 + p) 2 



+ 



4e 



1 + Pa 



2 + Pa 



U'a 



p a = (2 + V3) 1 / 3 + (2 



1 -5/? Q 

( 

Pa 

^) 1 / 3 - : 



(18) 
(19) 
(20) 
(21) 



(22) 



The combination Xp remains equal to 1 along the renormalization flow. In equa- 
tions (|TJ)-(22) all p-derivatives are at fixed bare parameters. The renormaliza- 
tion group flow has three nontrivial fixed points: the well-known directed per- 
colation fixed point with v = w = and u = nop = 2e/3, the symmetric 

(w = 0) fixed point (u s ,v si p s ) = (2e, 27e/64, 1/2) and the asymmetric fixed 
point (u a ,v a ,w a ,p a ) with 



(23) 



at leading order in e. The continuous phase transitions described by these fixed 
points have already been studied in other publications jl0| |U| ]. The symmetric 
fixed point (w = 0) is unstable with respect to the variable w. It corresponds to the 
case of equal diffusion constants Da = Db, whereas the asymmetric fixed point 
with w < governs the critical behavior for Da < Db- Since the sign of w is 
conserved along the renormalization group flow the asymmetric fixed point cannot 
be reached for w > 0. Therefore there is no fixed point for w > (Da > Db). In 
order to study the phase transition for Da > Db, which is the regime of interest, 
we consider in the next sections the solutions of the renormalization group flow in 
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more detail. Figure [j] shows a plot of the steady-state density of B's as a function 
of r (the deviation of the total density with respect to its mean-field critical value), 
for A = 1 and A > 1. 

Assuming the existence of a first-order transition for A < 1 (which would make 
of A = 1 a tricritical point), the jump of the order parameter across the transition 
point is related to the properties of the symmetric fixed point and should scale ac- 



cording to the tricritical scaling predictions developed by Lawrie and Sarbach [16], 



PB (r-)-p B (T+)<xay s , <* = -— ^- (24) 

d + js 

where 7 S is the Wilson function 7 evaluated at the symmetric fixed point. 



3 A phenomenological theory 

3.1 What happens when D B < Da? 

In dimension d < 4 the renormalization group flow has no stable fixed point at 
finite coupling constants. Nevertheless, we still expect a phase transition. Here 
follows a heuristic argument leading to the conclusion that this is a first order tran- 
sition. It is based, essentially, upon adding to the mean-field equations Eqs. <J3[> 
in an approximate way the fluctuations in the A particle density. Several steps in 
the argument are open to criticism but we expect it to provide the right qualitative 
picture. 

Let us consider the system at total particle density p and write 

P = Pc + Po (25) 

where p c is the critical density. The mean-field values of the stationary A and B 
densities are p\^ = p c and p% = po, respectively. 

We imagine the system divided into regions ("blocks") of volume L d , where L 
is arbitrary. Consider a particular block. The instantaneous density in this block is 
a fluctuating variable that we denote by 

PL = Pc + Po + $PL (26) 

where 5pi is a random term of average zero. 

We present the argument for the case po ^ Pc, ie- the average B density is 
much smaller than the average A density. Then the fluctuations of the total density 



are practically identical to those of p A - We have in particular (Sp 2 L ) = p c L d , so 
that the probability distribution of bpi, is 

P(^ L ) = Cexp(-^^) (27) 

A density fluctuation 5pi will relax to zero diffusively, hence on a time scale 

L 2 
7fl,L ~ — (28) 

We are now interested in fluctuations of pi well below the critical density ("nega- 
tive fluctuations"), say less than p c — p\. We have 

Prob(p L < p c - Pl ) ~ exp ( - L \ + pc Pl? ) (29) 



Such a fluctuation will still have the decay time Tfl £, given by (|28j) and therefore 
stay negative during a time 

T ne g ,L ~ ; -pr~, (30) 

In the meanwhile the local density of B particles will tend to zero with a relax- 
ation time T re i : £ which, according to the mean field equations, in the absence of B 
diffusion is given by 

T Td , L ~ y^- (31) 

kpi 
If pi is so large that T m \^ < T negi i, then during the lifetime of the negative density 



fluctuation the B particles will become locally extinct. Upon combining (30) and 
( pl| ) we find the condition for such a B extinguishing density fluctuation. We now 
take pi exactly large enough for this condition to be satisfied, but not larger, since 
we want to take into account all extinguishing fluctuations. This leads to a relation 
between p\ and L, viz. 

/i .1 (32) 



po + pi kL 2 

We now ask what the typical time interval T- mi l is between two such fluctuations 
in the same block. A rough estimate can be made as follows. The fraction /_ 
of time spent by the fluctuating density pi below the value p c — pi is equal to 



/_ = Prob(px < p c — pi), hence given by (g9|). This fraction is composed of short 

10 



intervals of typical length T neg) £ given by (|30D. The short intervals are separated by 
long ones of typical length T mi ^ that make up for the remaining fraction, 1 — /_, 
of time. Hence T negiL /T intiL = /_/(! - /_). Using © and © we find 



'int,L 






exp 



2p 



(33) 



where our replacing the prefactor pi/{po + pi) is without consequences for the re- 
mainder of the argument. The quantity T\ nt ^ is the decay time of the B population 
due to density fluctuations on scale L, in the absence of B diffusion. We now take 
into account the effect of this diffusion. The time Tb,l needed for a B particle to 
diffuse over a distance of order L is 



L 2 
TB > L ~ D~ B 



(34) 



All B particles will be eliminated from the system by negative density fluctuations 
on scale L unless Tb,l ^5 ^Int.L- By comparing (33) and (|34|) we obtain for the 
existence of a stationary state with a nonzero B density the condition 



t(T ,_ L d (p 1+Po ) 2 D A 

f( L , Po) = = > In -=- 

2p c D B 



for all L > a 



(35) 



where a is the lattice parameter and with p\ related to L by (32). The key point is 



now that when Db < Da, the inequality fl35[ ) can be satisfied only for po above 
some threshold poc to be determined below. Therefore 



P'c 



Pc + POc 



(36) 



is the new critical density. Since after having survived a negative density fluctuation 
any local B particle density will rapidly return to its average value p , there is at 
the new critical density a jump in p B equal to 



APS 



POc 



(37) 



We now determine the threshold value pq c . Since the inequality (35) has to 
hold for all L, we first determine the minimum value of its LHS as a function of 
L. In practice the calculation is most easily done by using p\ instead of L as the 
independent variable. The minimum occurs for L = ^ m ; n with 



£ 2 - 

smin 



(2d)- 2 (4 + d)(4-d)^- 
kp 



(38) 
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The values of p\ and f(L; po) at L = £ m j n are 

2d 

Pl.min = 1 iPO (39) 

4 — o 
f(U*Po) = (2 Pc )-iC d (A- d )-^(^y /2 p t-^ (40) 

where C d = (2d)~ d (4 + d) 2+d/2 . The condition £ min > a leads to 

(4-d)-^>l (41) 



and can always be satisfied by choosing a small enough. Upon inserting (40) in 
(|35|) we find the critical value po c below which there cannot exist a phase with B 
particles: 



AfD A \ 4- d/ _ _._/>, 



p 0c = est. pt a l-y) (4 - d) ln4-d -^ (42) 

Consistency requires that (^Tj) be satisfied when for po we substitute po c taken from 
(p2[). This leads again to a condition that can always be satisfied for a sufficiently 
small, whatever the dimension d. Let £* be the value of £ m ; n at which the existence 



condition Eq. ( p5[ ) of the B phase gets violated when po — > Pq c . One readily finds 



e,= (%) 4!d ln-^^ (43) 

This is the spatial scale at which the instability sets in that causes the first order 
transition. We also have to check that po c -C p c , in order to be consistent with po <C 



p c which was assumed following Eq. (26). This condition is certainly satisfied in 



the limit Dg — > D A , that we shall consider now. Setting as before 



o- = i - -=e (44) 



we obtain from the preceding equations 



2 

p 0c ~ est. (4 - d) a^-d (a -► 0+) (45) 



2 

e*-(^ ra a-ra (a-0+) (46) 



A:p ( 1/2 
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The relaxation time towards zero of the B density in the vicinity of the new critical 
density is 

Tb = T B ^ 

4 

a 4-d (cr^0 + ) (47) 



rv_/ 



Comparison of Eq. (45) with the tried tical scaling predictions of the previous sec- 



tion leads, with e = 4 — d, to the identification 

5 = \ (48) 

We expect that the exact theory gives power laws for the same quantities as the 
heuristic theory does, although with different exponents. One reason for this is the 
difficulty of correctly keeping track of the lattice parameter a. 

3.2 A nucleation picture 

Finally we would like to draw a parallel between the heuristic arguments developed 



above and the kinetics of first order transitions [ ]17| ] in thermodynamic systems. In 
those systems the description is based on a nucleation picture : the transition from 
a metastable to a stable phase occurs as the result of fluctuations in a homoge- 
neous medium. These fluctuations permit the formation of small quantities of a 
new phase, called nuclei. However the creation of an interface is an energetically 
unfavored process, so that below a certain size nuclei shrink and disappear. Nuclei 
having a size greater than a critical radius £* will survive an eventually expand. 
The analysis of the competition between bulk free energy and surface tension leads 
to an estimate of a critical nucleus size. 

In the original reaction-diffusion problem there is of course no such concept 
as a bulk free energy or surface tension. Nuclei are analogous to regions that are 
free of B particles. Those analogies should not be overinterpreted: they merely 
reinforce the intuitive picture of the reaction-diffusion processes leading to the 
first-order transition. 

4 Simulations in two dimensions: a hysteresis loop 

In this section we present the results of simulations performed on a two-dimensional 
500 x 500 lattice with periodic boundary conditions. At the beginning of the sim- 
ulation particles are placed randomly and independently on the sites of the lattice, 
with an average density p = 0.2. The ratio of the B particle density to the total 
density is arbitrarily chosen equal to 0.3. The decay probability of the B particles 
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is 7 = 0.1 per time step, and the contamination probability is k = 0.5 per time 
step. These parameters are held fixed. The diffusion constants Da and Db are 
varied. 

In each time step the reaction-diffusion rules are implemented by the following 
three operations. 

1 . Each B particle is turned into an A with probability 7. 

2. Each A (B) particle moves with a probability AD a (4Db); a moving parti- 
cle goes to a randomly chosen nearest neighbor site. 

3. An A particle is contaminated with probability k by each of the B particles 
on the same site. 

Then the new value of the average B density is evaluated and a new time step is 
begun. The process stops either when the system has fallen into its absorbing state 
or when the B density appears to have stabilized "(active state)". The latter sit- 
uation is considered to be reached when the slope of PB{t), as measured from a 
linear fit to the last 100 time steps, is 10~ 5 times as small as the maximal variation 
of the density of those 100 points. 

After this fixed density run we construct as follows a starting configuration for 
a new run in which the total density is increased by a factor 1.004. Two situations 
may occur. If at the end of the run just terminated the B particle density is nonzero, 
then we obtain the new starting configuration from the final one of the preceding 
run by randomly placing extra particles on the lattice while keeping the ratio of B's 
to the total number of particles constant. If at the end of the run just terminated the 
B particle density is zero, then we construct a new starting configuration with a B 
density equal to its value in the starting configuration of the preceding run. 

A new run is carried out at the new density. This process is iterated until some 
upper value of the total density is reached, taken equal to p = 0.5 in the present 
simulations. After that we carried out a step-by-step decrease of the total density, 
using a reduction factor of 0.996 per step, until we reached again the total density 
p = 0.2 of the beginning of the simulation. 

This whole procedure constitutes a simulation cycle. In this way we produced 
21 cycles with different pseudo-random numbers. Figures ^ and ^] show the result- 
ing order parameter curves in two different cases: in Figure || we have Db < Da 
and the system cycles counterclockwise through a hysteresis loop, signalling the 
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occurrence of a first order transition. In Figure [3] exactly the same simulation pro- 
cedure does lead to some dispersion in the order parameter curves, but not to a 
clearcut hysteresis loop; in this case the transition is known to be continuous. 

5 Perturbative calculation of the equation of state 

5.1 One-loop perturbation expansion 

In this section we determine the one-loop equation of state. We start from the 
dynamic functional of Eq. (0), in which we have included a particle source term 
— J d d rdt \h(r, i)V>(r, t). The equations of motion for the fields are 

= — = d t m - X(pV 2 m + aV 2 ib) + 2XpV 2 m (49) 

dm 

= -j = dtip + A(r - V 2 + fm)ip + ^ 2 - \gi/>$ - Xh (50) 

dip 2 

Note that the source term is not necessarily constant. Equations (G^) and (BO) are 
valid when they are inserted in averages. Taking the averages of ( |49[ ) and ( pQ ) we 
find for the densities M(r) = (m(r, t)) and ^(r) = (ip(r, t)) in a stationary state 
the exact equations 

P V 2 M(r) + (rV 2 $(r) = (51) 

T - v 2 + /M(r) + |*(r)] *(r) + fC m ^r) + |c^(r) = h(r) (52) 

with the correlation functions 

a™/,(r) = ((m - M(r)) (^ - ^(r))) C^t) = ((^ - ^(r)) 2 ) . (53) 

^From equation ( |5"Tj ) it follows that 

M(r) = --*(r) + c(r) (54) 

P 

where c(r) is a harmonic function (V 2 c = 0). Here we assume that c is constant. 
(It has to be constant in the thermodynamic limit if both M and \P are free of 
singularities and finite for r — ► oo.) If J v d d r m(r, t) = (which can always be 
achieved by a shift of r) we get 

c = ^7 / d d ^(r) (55) 

where V denotes the volume of the system. 
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The mean field equation for the profile reads 



r + /c-V 2 + f^ m/ (r) * m/ (r) = /i(r) 



(56) 



where g = g — 2a f /p. Equation ( pq ) shows that stability of the mean field theory 
requires that g > 0. For negative g higher powers in ifi have to be taken into 
account in the functional S. We first consider an external field h which is constant 
within a sphere of radius R and vanishes for r > R. For simplicity we take the 
thermodynamic limit V, R — > oo in such a way that R d /V — ► 0. In this case the 
region outside the sphere acts as a reservoir for the homogeneous mode of the field 



m, and in (56) we can set c = for r > 0. 



To calculate the one-loop correction to the equation of state we shift the fields 
tp and m by their average values and obtain 

S[fn,M + m;ip,ip + ^] = 5 [^;#] + S G [m,m;ip,ip;^] + Si[m; $, ijj\ (57) 

with 



S $;*] 



dt / d a r\ifj 



T- V 2 + -f) * 



(58) 



where we have kept the full r-dependence of VP and expressed M in terms of \t. 
The Gaussian part of S reads 



Salm, m; tp, ip; \P] = / dt / d d r [m (ftm - XV 2 (pm + at/))) - Ap(Vm) 5 



V> (ftV> + A (r - V 2 + (g - af/p)V) ij> + \f*m) - -\g^f z 
and the interaction part is 



S/fmjV',^] 



dt / d d r\ 



fmipip H — iftigip ~~ 9' l l>)' l l , 



(59) 



(60) 



The Gaussian propagators G,/, = (V^). G m = (mm), G m ^ = (imp) and G^, m = 
(i/jfh) follow from (^). They satisfy the differential equations 

(ft + A(f - V 2 )) G f (r, r'; i - t') + A/tfG^r, r'; t - t') = 6& (r - r')S(t - t') (61) 

(ft - XpV 2 ) G m ^(r,r';t-t') - AaV 2 G^(r,r';t - t') = 0(62) 

(ft - XpV 2 ) G m (r,r';t - t') - AaV 2 G^ m (r,r'; t - t') = 6^(r - r')S(t - i') (63) 

(ft + A(f - V 2 )) G^ m (r,r';t- t') + A/^G m (r,r';i- t') = 0(64) 
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with t = t + (g — af/p)ty. The Gaussian propagators can be used to determine 



the equal time correlation functions ( |53[ ) to lowest nontrivial order: 



*t /dV [A5^(r , )G m ^(r,r';t)G^(r,r';t) 
o J 

+2\p(VG m {r,r'-t)) (V'G^ m (r,r' ; t))] (65) 



C^(r) = / d£ / dV 



A^r') (G^(r, r'; t)) 2 + 2Ap (V'G^ m (r, r'; t)) 2 



.(66) 



For constant ^ the equations (|6T|)-(|64|) can easily be solved by Fourier transforma- 



tion. In this way one obtains for the fluctuation term in the equation of state (52) 



J^rm(! + n^i/> 



* 



f+(l+p)q 2 



1_ 2 g(gpg 2 + 2/ 2 fr) 

4 55 J + 4(r'+g 2 ) 



(67) 



with r' = r — af^/p. After dimensional regularization the momentum integral 
becomes 



J^mtp + T^V' 



A e tf 



l-e/2 



~gg - 4/ 2 + 55P + ^g(gpr' - 2f 2 *) 

T 



2e{l + p) \l + p 

e (l + p)V g(gpr'-2fH)^ (I + p)r 



+- 



2 r f - (1 + p )t' 

where A £ = (4^)~ d / 2 r(l + e/2)/(l - e/2). 



In- 



+ 0(e 2 ) 



(68) 



5.2 Renormalized equation of state 

In order to absorb the e-poles in the equation of state (|2,68) into a reparametriza- 
tion of coupling constants and fields we make use of the renormalized quantities 
tp R , A R , p R etc. introduced in Eqs. (11 Jl~2). 

The renormalized quantities satisfy the equation of state 



"> = **< t ' + T*' + WTk) 



— —2 - 

((1 + p R )u -Av- 2w) -^— In f-^ (69) 
1 + p R 1 + p R 



(p R u - 2w)t' r - 2vg R m R ( _ p 2 r R , _ 2 , 

• — — I r R In — (1 + p R )r R \n(p r R ) 



r R - (1 + p R )r R 



I + Pr 



where g R = g R - 2f R /p R , Tr = t r + (g R - f R /p R )^ R and r R = r R + c/ R ^ R are the 
renormalized counterparts of g, f and r', respectively. For / R = v = w = we 



recover the one loop-equation of state for directed percolation [20]. To simplify 
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the writing in equation (|^) the geometrical factor A e , the momentum scale p and 
a R have been absorbed into a rescaling of ^ R and h R , i.e. 

^V /VA - * R and ^ £ - 1 /V £/ V R / lR - /i R (70) 

Hereafter we will drop the index "R" since only renormalized quantities will be 
used. 



6 Flow equations 

6.1 Renormalization group for the equation of state 

The renormalizability of the field theory implies a set partial differential equations 
for the vertex functions. These are the renormalization group equations which 
follow from the independence of the bare vertex functions on the momentum scale 
p. To investigate the equation of state in the critical region we need the one-point 
vertex function jS 1 ' ' which is (up to a factor A) equal to h(r, VP; u, v, w, p; p). The 
renormalization group equation for h reads 



d d d d d d 

<-TT- + PuT- + PvTT + Pw^~ + PpTT + «T— - +C~7 K T 1 *5 U i V i W ■> P'l MJ = ° 
Op OU OV OW Op OT 

(71) 
with the /3-functions given in Sec. 2. 

6.2 Scaling form of the equation of state 



The renormalization group equation Eq. < J7 1| ) can be solved by characteristics with 
the result 

h(r, *; u, v, w, p; p) = Y h (£)-\p£) 2+d l 2 
xh (Y T (e)(p£)- 2 T, (pl)- d / 2 *; u(£),v(£), w(£),p(£); l) (72) 



where 

do(I) 



din 

dlnY h (£) 



P a (u(£),v(£),w(£),p(£)) for a = u,v,w,p (73) 

K (u(£),v(£),w(£),p(£)) (74) 



dln-l 

dlnY T (£) 



din j 



7 (u(£),v(£),w(£),p(£)) - C(u(£),v(£),w(£),p(£)) (75) 
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are the characteristics with the initial conditions u(l) = u, v(l) = v, w(l) = w, 
p(l) = p and Y T {\) = l/i(l) = 1. Some solutions of the flow equations ( f73| ) are 
depicted in figure f|. For small initial values of w a Da — Db the trajectories 
first approach the unstable manifold of the symmetric fixed point (w = 0) before 
they flow away. The unstable manifold leaves the stability region of the mean field 
theory (g > 0, or u — 2w/p > 0) at the point (u*, u*, u;*,/?*). The numerical 
solution of the flow equations yields u* = 12.32e, v+ = 2.104e, w+ = 22. lie, 
and p^ = 3.589. The intersection point of the unstable manifold with the stability 
edge g = is of special interest since the perturbatively improved mean field 
theory should be a good approximation for small g (as will be discussed in the next 
section). Therefore the phase transition for very small w > is governed by the 
equation of state at the point (u*, v*, W+, p*). In the following we shall denote by Z± 
the value of the flow parameter at which «(£*) — 2w{i^) j p{i^) = and by £* = e e * 
the related length scale. One can identify £* with its heuristic counterpart defined 



in Eq. <]43t)- The flow equations are too complicated to be solved analytically for 
all £, however it is possible, using scaling arguments, to predict that, as a —> + , 
£* oc <t 2 / 7s , a form also proposed in Eq. ([46]). 



7 A first-order transition for Da > D B 

In this section we study the mean field equation of state with the one-loop fluctua- 
tion correction (^9)) for small values of the coupling g -C g- Our motivation is an 
analogy of our reaction diffusion system with spin systems with cubic anisotropy. 
Using a modified Ginzburg criterion for systems with cubic anisotropy Rudnick ^ 
(see also [^]) has shown that the fluctuation corrected mean field approximation 
should give reliable results if the values of the coupling coefficients are close to the 
stability edge of the mean field theory. In the case of the reaction diffusion system 
with Da > Db the stability edge is given by g = 0. 

Now assume that we start with a very small value of the coupling w and choose 



the flow parameter £ = £* in (72). After the ^-dependent prefactors Y T (£i,)£~ 2 
etc. have been asborbed into a rescaling of r, h, and ^ the improved mean field 
equation of state takes the simple form 



h = V 



Ui,-4Vi, / g* m \, r + g^ 2 

T + J7TT — ^ T + V^ ln 2/i , — r + O two-loop 
4(1 + pi,) 2 V 2 / p z (l + pi,) 



(76) 



with (u* - 4w*)/(4(l + p^) 2 ) = 0.04635 e > 0. 

In the limit h —> + the absorbing state with ^> = is a solution of the equation 
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of state for all r. For r < T sp i no d with 

t sp w = /u e — — r = 0{e) (77) 

4(1 +p*) 

there is also a solution with ^ > and d^> jdh > (see figure |5|). To see this 
one should anticipate that when the order parameter is of the order of its value at 
r S pinod one has 

Tspuiod = 0(e) (78) 

and to leading order the reasoning is carried out on 



fc = ¥ 



u*-4t>* g* g^ 

T + -T7Z ^r-r* In 



4(1 +p*) 2 2 2(l+^)/x 2 



(79) 



Solving the latter equation for h = and ^4 = yields, to leading order in e, 



_ 8(1 + p±) 2 Tspjnod _ 2 -1 U* ~ 4v* 

* spinod , ? ^"spmod [A & A , , l°Oj 

u* - 4^* ^ ^ 4(1 + p*) 



which a posteriori justifies working on the approximation Eq. (|79j). Since the sus- 
ceptibility x = d^/dh (the response of the order parameter to particule injection) 
is positive the new solution is at least metastable for all r < T sp i no d- The open 
question is for which range of r the solution is stable, i.e., stable also with respect 
to nucleation processes. If we could derive the equation of state from a free energy 
this question would be easy to answer: a metastable solution is a local minimum of 
the free energy while a stable solution is a global minimum. However, since there 
is no free energy for the reaction diffusion system we have to look for a different 
stability criterion, e.g. the form of density profiles. By analogy with equilibrium 
systems we expect that the solution with ^ > is stable below a coexistence 
"temperature" (in our case we should talk about a coexistence density) T coex with 

U ^- T coex < Tspinod- 

8 Deriving r coex from a density profile 

If r is so large that equation ( f76| ) has only the trivial solution ^ = the density 
^(r) generated by a local source h(r) decays rapidly with increasing distance from 
the region in which h(r) is nonzero. Consider for instance a plane particle source 
with 

h(r) = h 5(r ± ) (81) 
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where r± is the coordinate perpendicular to the source. For simplicity assume 
that ho —> oo. For large r± the profile [i.e., the solution of (52)] will become 
independent of r± and tend to either (a) ^ = if r is sufficiently large or (b) the 
nonzero solution of (f76|). In the thermodynamic limit the whole profile \P(r) is 
uniquely determined by the source h(r) since ^> has to be finite for r± — ► ±00. In 
the case (a) we can conclude that r > T coex whereas (b) occurs for r < T coex . 

In order to calculate a density profile perturbatively one usually starts with the 
mean field profile and assumes that the fluctuation corrections are small (of the 
order e, say). In the present case this procedure would not lead to the desired 
result since the nonzero solution of the equation of state for r < T spino( i is of the 
order e . It can therefore not be derived as a small correction to ^ m f- Instead we 
have to compute the equation for the density profile perturbatively and then study 
the asymptotic behavior of its solutions. This means that we need the correlation 
functions C m ^,(r; {*£}) and C^(r;{^}) in equation (52) for a general function 
\f(r). Of course, we cannot compute these functions exactly but it is possible to 
derive a systematic gradient expansion for C mi p and Cw,. 



We first compare the leading e-orders in equation ( |52| ) (with M = —aSSf/p): 
Since we expect that the coexistence point to be located in the parameter range 
< r < T spin0( i = 0(e) we may set r = 0(e). The limit of a very weak first 



order transition is governed by the point 
and that the combination fC m j, + 



u±,v+,w+ 



, p*) which means that g = 



3.(1 , 

72v 



is of order 0(e). Comparing the terms on 
the l.h.s. of (|2|) therefore yields (V 2 ^)/^ = 0(e), i.e. gradients of * may be 
considered as small quantities when we calculate C m ^ and Cw,. Using the Taylor 
series 



Hf(r')=^(r) + J2^ £ (r'-rU...(r'-r) aN d ai ...d aN y(r) (82) 

N=l ' ai—ctN 

for the profile one arrives at a gradient expansion for C m ^ and C^ of the form 

oo _. 

C^(r,{¥}) = C^(¥(r)) + £— £ C^ ai ... aN (V(r))d ai . ..0„„¥(r). 



N=l ai—acff 



(83) 



At leading order in e only the first term on the r.h.s. of (83) contributes, i.e., we may 
simply replace ^ in (^) by the profile ^(r) and use the result in equation (52). 
After application of the renormalization group as before and absorbing ^-dependent 
prefactors the equation for the profile becomes 



h(r) = *(r) 



r + 



4(1 + P*) 



T + ^^( r ))l n 



r + g^(r)/2 
/^ 2 (1 + P*) 



V 2 ^(r) + 0(e 2 ) 
(84) 



21 



In order to extend this result to the next order in e one has to (i) compute C m ^ and 
Cw, for constant ^ to two-loop order and (ii) take into account the V 2 ^ -correction 



in the gradient expansion (Q) to one-loop order. In this way the V^-term in (34) 
may receive a ^-dependent correction. 

For h(x) = h 5(r±) Eq. (Q) can be integrated once after multiplication of both 
sides with i ff'(r±). In figure gthe result is depicted for various values of r. There 
is a value T coex < T sp i no d such that the profile does not tend to zero for r± — ► ±00 
if r < T coex . As discussed above T coex is the coexistence point below which the 
active phase becomes a stable solution of the equation of state. Again one works 



with Eq. (79). In practice one has to solve the system composed of Eq. (79) with 
h = along with its integrated counterpart 







T, 



coex T 2 



» + 



Av. 



2 coex 24(1 + p*) 2 



9**\ 



In 



g^c 



1 



2m 2 (1 + P* 



(85) 



The explicit calculation yields ^ coe x = L izf^ \ coex (which also justifies the use 
of the simplified equation of state Eq. (79)), with 



T _ „2 -2/3^*. 

'coex — f* c 



4u* 



u.yoT S pj no ^. 



6(1 + p 

Tcoex /T sp i n od is not a universal number, but the susceptibility ratio \+/X- with 

a* 

X_ = lim 



(86) 



Xh 



km — - 

h-^Q Oh 



X- = lim — — 
/i— 5/i 



(87) 



*>o, 



is universal. To one-loop order one finds 

X- 2 



+ C?(e) 



(88) 



This result is analogous to the universality of the magnetic susceptibility ratio 
found by Rudnick [Q] and Arnold and Yaffe []3]]. 

9 Conclusions and prospects 

9.1 A heuristic functional for the steady state phase diagram 

In this paragraph we would like to build a posteriori a functional of the order 
parameter field \t (r) describing the phase diagram in the stationary state of the 
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system. We emphasize that the following is only valid to one loop order. We define 

F{^] by 



(89) 



F[V] = fd d r i(W f + f dV r (1 ' 0) [0, -0] 

By construction of course |£ = is equivalent to the equation of state Eq. 
It is instructive to plot F as a function of vl/ for various values of r (see figure [7]). 
There it appears possible to deduce the phase diagram from the global minima of 
F\$>]. However the route leading to the functional F[\P] follows a series of field- 
theoretic detours. The suggestive notation F, which reminds of a free energy (in 
the equilibrium statistical mechanics sense) is however misleading. For instance 
it could not be used as an effective Landau hamiltonian for the calculation of a 
Gibbs partition function describing fluctuations directly in the steady state. This 
functional is a remarkably compact and intuitive way of summarizing the properties 
of the steady state phase diagram. In particular the spinodal point T sp i no d appears as 
the point below which F develops a second minimum. Below T coex that minimum 
becomes the global minimum. To one loop order the equilibrium vocabulary can 
therefore be used carelessly. 

9.2 Summary 

In the course of this work we have elaborated a phenomenological description of 
a fluctuation-induced first-order transition taking place in a nonequilibrium steady 
state, the first one of this sort. We have in parallel applied field-theoretic tech- 
niques to derive an effective (renormalization-group improved) equation of state 
that incorporates those fluctuations. This yields a PDE for the order parameter 
in the steady state. Performing a study of the stability (against space fluctuations 
of the order parameter) of the solutions of this PDE has led us to a complete de- 
scription of the phase diagram. We have identified in this nonequilibrium situation 
a concept analogous to the point of spinodal decomposition consistent with our 
phenomenological description. 

9.3 Possible applications 

Among the many nonequilibrium systems that appear in the literature, driven dif- 
fusive systems lend themselves to an analytic treatment by techniques similar to 



those of the present article [21]. In a number of such systems though, a significant 
portion of the phase space (in terms of control parameters) escapes conventional 
analysis. In some cases we believe that the reason is the occurrence of fluctuation- 
induced first order transition, such as in [^] or [g|. It would be quite interesting to 
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see how both the technical argument and the heuristics can be extended to those 
systems. This will be the subject of future work. 
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Appendix 

There are two ways to proceed in order to obtain the equation of state to one-loop 
order. In this appendix we follow the route familiar from static critical phenomena. 
We use the action Eq. (||) as the starting point of our analysis. The first task is to 
determine the one-loop expression of the effective potential T. 

1 C d^o duo 
T[^]=S[i/i,if>] + -J —jL — ]n6etS ,, (q,u)$,il>] (90) 

with the matrix S" defined by 

/ 5 2 S 8 2 S \ 



S"(q,u)=\ W-1'^W.") <5<K-q,-w)^(q,w) (91) 

\ 5i/)(— q,-aj)5i/)(q,Lij) Sip(— q,— ui)5tp(q,ui) / 

S'diq, u) = A(q 2 +f)-iu- \g* - \ 2 af*^— - 2^4^ (92) 

q z — iuj q 4 + u> z 

SU<l,") = -^-4* 2 -W (93) 
q 4 + to z 

S»(q, w ) = -Hj^^-2*^ + A^ (94) 
q 4 + or g 4 + uj z 

S£ 2 (q, u) = A(q 2 + f) + iu; - XgV - A V*^— - 2%^^ (95) 

q z + zw q 4 + ur 

For a homogeneous source term h the equation of state for a homogeneous order 
parameter ^ now follows from the requirement that 

-£[0,*] = (96) 

dtp 

It is a tedious but straighforward task to find the one-loop correction to T^ 1,0 ) in the 
form of an integral over momentum and frequency. Writing 



r (1 ' 0) [i/> = o, ip = *] = -\h + \^{t + -g&) + jr (1 ' 0) 



(97) 



one finds 



x(5(q 4 + W 2 ) + 2A/ 2 q 2 *) (98) 



x [(w 2 - iAw - B) (a; 2 + »Ai> - B)] 
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where we have defined the auxiliary variables 

,4 = (A + l)q 2 + Af , B = q 2 A(q 2 + t') 
Upon using the following integration formulas, 



duj 



1 



2vr uj 2 ± iAuo - B 
Eq. ( |98| ) simplifies into 

^(1,0) = * f d d q 



0. 



duj 



1 



2ir \uj 2 ± iAcv - B\ 2 2AB 



(99) 



(100) 



2 7 (27r) d 



^ +f A 9 * ) {; 9 ^ + T^ ) + Tr P * 



* + w +I ^> 



(101) 



The above expression can be cast in a form suitable to perform the q-integrals: 

d d q 



8T^)=-\xgj 



+ 



+ 



(27r) d 
A^ 



4((1 +p)r'-f) 

A* 

4(l + p)((l+p)r'-f) 

- 4p/ 2 f + 2g(l + p)/ 2 * 



[£(p<7 - 2cr/)r' - 2/ 2 s*] / 



d d q 1 



(2vr) rf q 2 + t' 



2gaff - (1 + p)gg^-V 



d d g 1 



(27r)^q2 + T ^ 



(102) 



We use dimensional regularization to compute the momentum integrals: 

d d q 1 2 n _ E/2 A f d d q 1 2 / f N 1 " t 2 



-V 1 " 6 / 2 ^ 



(2vr) d q 2 + tJ- " eVl + p' 



In terms of renomalized quantities the equation of state has the form 



(103) 



/'«■ %U R + ^^ R + ' 



+ 



2 4(1 + p R ) 

(p R u - 2w)t r - 2vg R V R 
f R - (1 + p R )r R 



((1 + p R )u - Av - 2w) -J- l n £j!l 
1 + p R 1 + p R 



r R ln 



p 2 r R 
1 + Pr 



^l + p R )r R ln(/i t r 



(104) 



which is Eq. (69) 



26 



Figures 



27 



Mt = +00)) 




D A < D B 



(3 = 1 




0(e 2 ) 



D A = £>f 



£>a > -Dje 



Figure 1: Phase diagram in the (T,ip(t = 00)) plane (the ordinate is the steady- 
state density of B particles) for A > 1, A = 1 and a conjecture for A < 1. Also 
shown is the order parameter exponent j3. 
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Figure 2: The order parameter ^ for AD a = 0.8, and ADb = 0.1, with 7 = 0.1 
and k = 0.5. The system cycles anticlockwise through a hysteresis loop. 
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Figure 3: Same as Fig. 2, but for AD a = 0.1, and ADb = 0.8. In this case the 
transition is known to be continuous and no hysteresis loop appears. 
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D A <D E 



D A > D e 



Figure 4: Flow diagram in the (gg, w) plane. The leftmost black dot stands for 
the Da < Db fixed point while the one lying on the w = axis stands for the 
symmetric Da = Db fixed point. They both describe second order transitions. 
Typical trajectories have been drawn. Those starting close to the symmetric fixed 
point but with an initial positive w eventually flow away from the symmetric fixed 
point. 
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7~ ^ 'spinod 



^ V 



Figure 5: Sketch of the function h(r, *), Eq. @. 
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^tf(rj_) 



*'(r ± ) 



'spinod 



Figure 6: Derivative of the density profile i &'(r±) as a function of VP(rj_) for r > 
Tcoex- The boundary condition at r± = + is given by ^'(0) = —ho/2. 
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T ^ ^spinod 
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spinod ^ T ^ Tcoe 




T *C T C oe. 



Figure 7: F\^i\ as a function of * for decreasing values of r. We identify the 
spinodal point as the value of r below which F develops a local nonzero minimum 
and the coexistence point as the point at which the two minima become degenerate. 
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